##### ################################################### ######
#####                                                     ######
#####   Input: clean experiment data                      ######
#####   Output: Benjamini-Hochberg corrected p-values     ######
#####                                                     ######
##### ################################################### ######

setwd("/Users/lotte/Dropbox/PhD/style_experiment/replication")
rm(list=ls())

# Load libraries

library(data.table) # CRAN v1.14.2
library(plyr) # CRAN v1.8.6     
library(dplyr) # CRAN v1.0.9 
library(tidyverse) # CRAN v1.3.1
library(ggplot2) # CRAN v3.3.6 
library(broom) # CRAN v0.7.11
library(patchwork) # CRAN v1.1.1
library(estimatr) # CRAN 0.30.6

# Load clean experiment data 

load("data/style_data.Rdata")
load("data/emotion_data.Rdata")
load("data/aggression_data.Rdata")
load("data/evidence_data.Rdata")

# Set seed for reproducibility
set.seed(1234)

# Run relevant models
emotion_perception <- lm(perceived_emotion ~ objective_style, data = emotion)
aggression_perception <- lm(perceived_aggression ~ objective_style, data = aggression)
evidence_perception <- lm(perceived_evidence ~ objective_style, data = evidence)
emotion_likeability <- lm(perceived_likeability ~ objective_style, data = emotion)
aggression_likeability <- lm(perceived_likeability ~ objective_style, data = aggression)
evidence_likeability <- lm(perceived_likeability ~ objective_style, data = evidence)
emotion_competence <- lm(perceived_competence ~ objective_style, data = emotion)
aggression_competence <- lm(perceived_competence ~ objective_style, data = aggression)
evidence_competence <- lm(perceived_competence ~ objective_style, data = evidence)
emotion_likeability_covariates <- lm(perceived_likeability ~ mp_gender*female_stereotype_congruent +
                                       respondent_gender + age + degree_educated, data = emotion) 
aggression_likeability_covariates <- lm(perceived_likeability ~ mp_gender*female_stereotype_congruent +
                                          respondent_gender + age + degree_educated,
                                        data = aggression)
evidence_likeability_covariates <- lm(perceived_likeability ~ mp_gender*female_stereotype_congruent + 
                                        respondent_gender + age + degree_educated, 
                                      data = evidence)
emotion_competence_covariates <- lm(perceived_competence ~ mp_gender*female_stereotype_congruent + 
                                      respondent_gender + age + degree_educated + political_attention,
                                    data = emotion)
aggression_competence_covariates <- lm(perceived_competence ~ mp_gender*female_stereotype_congruent + 
                                         respondent_gender + age + degree_educated + political_attention,
                                       data = aggression)
evidence_competence_covariates <- lm(perceived_competence ~ mp_gender*female_stereotype_congruent + 
                                       respondent_gender + age + degree_educated + political_attention,
                                     data = evidence)
emotion_prevalence_covariates <- lm(perceived_emotion ~ mp_gender*objective_style_emotion +
                                      respondent_gender + left_right_placement + age,
                                    data = emotion)
aggression_prevalence_covariates <- lm(perceived_aggression ~ mp_gender*objective_style_aggression + 
                                         respondent_gender + left_right_placement + age,
                                       data = aggression)
evidence_prevalence_covariates <- lm(perceived_evidence ~ mp_gender*prevalence_evidence + 
                                       respondent_gender + age + degree_educated,
                                     data = evidence)
  
p_values_all <- c(summary(emotion_perception)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(aggression_perception)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(evidence_perception)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(emotion_likeability)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(aggression_likeability)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(evidence_likeability)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(emotion_competence)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(aggression_competence)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(evidence_competence)$coefficients["objective_styleTreatment","Pr(>|t|)"], 
              summary(emotion_likeability_covariates)$coefficients["mp_genderWoman:female_stereotype_congruentCongruent", "Pr(>|t|)"], 
              summary(aggression_likeability_covariates)$coefficients["mp_genderWoman:female_stereotype_congruentCongruent", "Pr(>|t|)"], 
              summary(evidence_likeability_covariates)$coefficients["mp_genderWoman:female_stereotype_congruentCongruent", "Pr(>|t|)"], 
              summary(emotion_competence_covariates)$coefficients["mp_genderWoman:female_stereotype_congruentCongruent", "Pr(>|t|)"], 
              summary(aggression_competence_covariates)$coefficients["mp_genderWoman:female_stereotype_congruentCongruent", "Pr(>|t|)"], 
              summary(evidence_competence_covariates)$coefficients["mp_genderWoman:female_stereotype_congruentCongruent", "Pr(>|t|)"], 
              summary(emotion_prevalence_covariates)$coefficients["mp_genderWoman:objective_style_emotionEmotional", "Pr(>|t|)"], 
              summary(aggression_prevalence_covariates)$coefficients["mp_genderWoman:objective_style_aggressionAggressive", "Pr(>|t|)"], 
              summary(evidence_prevalence_covariates)$coefficients["mp_genderWoman:prevalence_evidenceAnecdote", "Pr(>|t|)"])

# Choose alpha level
alpha <- 0.05

# Without any corrections
sig <- p_values_all < alpha

# Conduct Benjamini-Hochberg procedure 
BH_sig <- p.adjust(p_values_all, "BH") < alpha

# Collect model names
model_names <- c("Unconditional emotion perception", "Unconditional aggression perception", "Unconditional evidence perception", 
                 "Unconditional emotion likeability", "Unconditional aggression likeability", "Unconditional evidence likeability", 
                 "Unconditional emotion competence", "Unconditional aggression competence", "Unconditional evidence competence",
                 "Conditional emotion likeability", "Conditional aggression likeability", "Conditional evidence likeability", 
                 "Conditional emotion competence", "Conditional aggression competence", "Conditional evidence competence",
                 "Conditional emotion perception", "Conditional aggression perception", "Conditional evidence perception")

# Add names, p-values, and significance into dataframe
adjusted_unadjusted <- data.frame(model_names,  p_values_all,  sig <- p_values_all < alpha, p.adjust(p_values_all, "BH"), 
                                  p.adjust(p_values_all, "BH") < alpha)

# Rename dataframe 
names(adjusted_unadjusted) <- c("model_names", "uncorrected_pvalues", "uncorrected_significant", "corrected_pvalues", 
                                "corrected_significant")

# Save data 
save(adjusted_unadjusted, file = "data/corrected_pvalues.Rdata")
